In [1]:
import plotly
import plotly.graph_objs as go
import plotly.express as px
from plotly.subplots import make_subplots

import torch
import numpy as np
import pandas as pd

import dill as pickle

import seaborn as sns
import random
import math
In [2]:
def get_rms(records):
    """
         Среднеквадратичное значение отражает эффективное значение, а не среднее
    """
    return math.sqrt(sum([x ** 2 for x in records]) / len(records))
In [3]:
grid_res = 70
title = 'wave_equation'
df = pd.read_csv(f'data/{title}/wave_sln_{grid_res}.csv', header=None)
u = df.values
u = np.transpose(u) # x1 - t (axis Y), x2 - x (axis X)

sigma_NR = 0.01
n = []

for i in range(u.shape[0]):
    n.append(np.random.normal(0, sigma_NR*get_rms(u[i,:]), u.shape[1]))
#     print(f' rms = {get_rms(u[:,i])}, std = {np.std(u[:,i])**2}, Standard deviation = {sigma_NR*get_rms(u[:,i])}')
n = np.array(n)
print(f' rms_all = {get_rms(u.reshape(-1))}, std_all = {np.std(u)**2}, Standard deviation = {sigma_NR*get_rms(u.reshape(-1))}')

u_total = u + n


# noise in other form
u_noise_all = []
variance = [0.01]#[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]

for v in variance:
    noise = np.random.normal(0, v, size = [u.shape[0], u.shape[1]])
    u_noise_all.append(u + noise)

# mesh    
x_grid = np.linspace(0, 1, grid_res + 1)
t_grid = np.linspace(0, 1, grid_res + 1)

x = torch.from_numpy(x_grid)
t = torch.from_numpy(t_grid)

grid = torch.cartesian_prod(t, x).float()

 rms_all = 3.868013079116736, std_all = 3.7768482896998234, Standard deviation = 0.03868013079116736
In [4]:
# building 3-dimensional graph
fig = go.Figure(data=[
    go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=torch.from_numpy(u).reshape(-1), name='Initial field',
              legendgroup='i', showlegend=True, color='rgb(139,224,164)',
               opacity=0.5),
    go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=torch.from_numpy(u_total).reshape(-1), name='Initial field + noise',
              legendgroup='i_n', showlegend=True, color='rgb(139,224,80)',
              opacity=0.8),
#     go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=torch.from_numpy(u_noise_all[0]).reshape(-1), name='Initial field + noise2',
#               legendgroup='i_n_3', showlegend=True, color='rgb(139,224,164)',
#               opacity=0.5)
])
fig.update_layout(scene_aspectmode='auto')
fig.update_layout(showlegend=True,
                  scene=dict(
                      xaxis_title='x1 - t',
                      yaxis_title='x2 - x',
                      zaxis_title='u',
                      zaxis=dict(nticks=10, dtick=1),
                      aspectratio={"x": 1, "y": 1, "z": 1}
                  ),
                  height=800, width=800
                  )
fig.show()
In [5]:
u_main = torch.load('data/confidence_region/file_u_main_35.pt')
prepared_grid_main = torch.load('data/confidence_region/file_prepared_grid_main.pt')
# u_main = torch.load('data/confidence_region/file_u_main_35_var_0.0001.pt')
# prepared_grid_main = torch.load('data/confidence_region/file_prepared_grid_main_var_0.0001.pt')
In [6]:
mean_arr = np.zeros((u_main.shape[1], u_main.shape[2]))
var_arr = np.zeros((u_main.shape[1], u_main.shape[2]))
s_g_arr = np.zeros((u_main.shape[1], u_main.shape[2])) # population standard deviation of data.
s_arr = np.zeros((u_main.shape[1], u_main.shape[2])) # sample standard deviation of data

for i in range(u_main.shape[1]):
    for j in range(u_main.shape[2]):

        mean_arr[i, j] = np.mean(u_main[:, i, j])
        var_arr[i, j] = np.var(u_main[:, i, j])
        m = np.mean(u_main[:, i, j])
        s_arr[i, j] = math.sqrt(np.sum(list(map(lambda x: (x - m)**2, u_main[:, i, j])))/(len(u_main[:, i, j]) - 1))

mean_tens = torch.from_numpy(mean_arr)
var_tens = torch.from_numpy(var_arr)
s_g_arr = torch.from_numpy(var_arr) ** (1/2)
s_arr = torch.from_numpy(s_arr)

# Confidence region for the mean
upper_bound = mean_tens + 1.96 * s_arr / math.sqrt(len(u_main))
lower_bound = mean_tens - 1.96 * s_arr / math.sqrt(len(u_main))

mean_tens = mean_tens.reshape(-1)
upper_bound = upper_bound.reshape(-1)
lower_bound = lower_bound.reshape(-1)
In [7]:
# building 3-dimensional graph
fig = go.Figure(data=[
    go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=torch.from_numpy(u).reshape(-1), name='Initial field',
              legendgroup='i', showlegend=True, color='rgb(139,224,164)',
              opacity=0.5),
    go.Mesh3d(x=prepared_grid_main[:, 1], y=prepared_grid_main[:, 0], z=mean_tens, name='Solution field',
              legendgroup='s', showlegend=True, color='lightpink',
              opacity=1),
    go.Mesh3d(x=prepared_grid_main[:, 1], y=prepared_grid_main[:, 0], z=upper_bound, name='Confidence region',
              legendgroup='c', showlegend=True, color='blue',
              opacity=0.20),
    go.Mesh3d(x=prepared_grid_main[:, 1], y=prepared_grid_main[:, 0], z=lower_bound, name='Confidence region',
              legendgroup='c', color='blue', opacity=0.20)

])

fig.update_layout(scene_aspectmode='auto')
fig.update_layout(showlegend=True,
                  scene=dict(
                      xaxis_title='x1 - t',
                      yaxis_title='x2 - x',
                      zaxis_title='u',
                      zaxis=dict(nticks=10, dtick=1),
                      aspectratio={"x": 1, "y": 1, "z": 1}
                  ),
                  height=800, width=800
                  )
In [9]:
fig.add_trace(go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=torch.from_numpy(u_total).reshape(-1), name='Initial field + noise',
              legendgroup='i_n_2', showlegend=True, color='rgb(139,224,80)',
              opacity=0.5))    
fig.show()
In [10]:
fig = go.Figure(data=
                go.Contour(x=prepared_grid_main[:, 0],
                           y=prepared_grid_main[:, 1],
                           z=mean_tens,
                           contours_coloring='heatmap'))
fig.update_layout(
    title_text='Visualization of the equation solution'
)
fig.show()
In [11]:
# title="Visualization of the variance"
fig = go.Figure(data=
                go.Contour(x=prepared_grid_main[:, 0],
                           y=prepared_grid_main[:, 1],
                           z=s_arr.reshape(-1),
                           contours_coloring='heatmap'))

fig.update_layout(
    title_text='Visualization of the variance'
)
fig.show()
In [12]:
import seaborn as sns
import random
import math
import statistics
In [13]:
with open(f'data/wave_equation/data_equations_1000_3.pickle', 'rb') as f:
    equations = pickle.load(f)
In [14]:
equations[-5:]
Out[14]:
[{'d^2u/dx2^2{power: 1.0}': 0.03969843028379964,
  'u{power: 1.0}': 0.23222484867099805,
  'C': -0.8529317724716573,
  'du/dx1{power: 1.0}': -0.006340683146567955,
  'd^2u/dx1^2{power: 1.0}_r': '-1'},
 {'d^2u/dx2^2{power: 1.0}': 0.05227599127289423,
  'u{power: 1.0}': 0.5139793960777936,
  'C': -2.1702427564824607,
  'du/dx1{power: 1.0}': -0.00888923877672482,
  'd^2u/dx1^2{power: 1.0}_r': '-1'},
 {'d^2u/dx2^2{power: 1.0}': 0.030759351423659955,
  'u{power: 1.0}': -0.24841363495537536,
  'C': -0.21197063543632294,
  'du/dx1{power: 1.0}': -0.012856358319632994,
  'd^2u/dx1^2{power: 1.0}_r': '-1'},
 {'d^2u/dx2^2{power: 1.0}': 0.023149952184001326,
  'u{power: 1.0}': -0.27905269468301036,
  'C': 0.799733208174493,
  'du/dx1{power: 1.0}': 0.006162742346482619,
  'd^2u/dx1^2{power: 1.0}_r': '-1'},
 {'d^2u/dx2^2{power: 1.0}': 0.03172777819846863,
  'u{power: 1.0}': -0.22339757514847,
  'C': -0.35161106435265055,
  'du/dx1{power: 1.0}': -0.013821403219482449,
  'du/dx1{power: 1.0}_r': '-1'}]
In [15]:
d_main = {}
k = 0
In [16]:
for i in range(len(equations)):
    for temp, coeff in equations[i].items():
        if temp in d_main:
            d_main[temp] += [0 for i in range(k - len(d_main[temp]))] + [coeff]
        else:
            d_main[temp] = [0 for i in range(k)] + [coeff]
        
    k += 1

for key, value in d_main.items():
    if len(value) < k:
        d_main[key] = d_main[key] + [0 for i in range(k - len(d_main[key]))]
In [17]:
d = pd.DataFrame(d_main)
d.columns[:]
Out[17]:
Index(['d^2u/dx2^2{power: 1.0}', 'u{power: 1.0}', 'C', 'du/dx1{power: 1.0}',
       'd^2u/dx1^2{power: 1.0}_r', 'du/dx1{power: 1.0}_r'],
      dtype='object')
In [18]:
for col in d.columns:
    d = d.astype({col: np.float64})
In [19]:
df_Nan = d#d.replace(0, np.NaN)
df_Nan.hist(column = d.columns[:], figsize=(20, 15), bins = 100, rwidth = 0.6)
None
In [20]:
df_Nan = d.replace(0, np.NaN)
df_Nan.hist(column = d.columns[:], figsize=(20, 15), bins = 100, rwidth = 0.6)
None
In [21]:
C_temp = d['d^2u/dx2^2{power: 1.0}']
C_2 = C_temp[C_temp != 0]
In [22]:
m = np.mean(C_2)
s = math.sqrt(np.sum(list(map(lambda x: (x - m)**2, C_2)))/(len(C_2) - 1)) # cтандартное отклонение выборки)
print(s)

0.011554727664242445
In [23]:
# Доверительный интервал для среднего при неизвестной дисперсии
maximum_estimation_error = 1.96 * s / math.sqrt(len(C_2))
print(f'maximum_estimation_error = {maximum_estimation_error}')
lower = np.mean(C_2) - maximum_estimation_error
upper = np.mean(C_2) + maximum_estimation_error
print(f'interval estimation: {lower} <= m <= {upper}')

maximum_estimation_error = 0.000717246117665491
interval estimation: 0.030428234878135833 <= m <= 0.03186272711346681